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Abstract 

o 

*^ ' General models of Gibbs Delaunay-Voronoi tessellations, which can be viewed as extensions of 

K^ . Ord's process, are considered. The interaction may occur on each cell of the tessellation and be- 

jrt ' tween neighbour cells. The tessellation may also be subjected to a geometric hardcore interaction, 

forcing the cells not to be too large, too small, or too flat. This setting, natural for applications, 
introduces some theoretical difficulties since the interaction is not necessarily hereditary. Mathe- 
matical results available for studying these models are reviewed and further outcomes are provided. 

C^ ' They concern the existence, the simulation and the estimation of such tessellations. Based on these 

results, tools to handle these objects in practice are presented: how to simulate them, estimate their 

Li^ ' parameters and validate the fitted model. Some examples of simulated tessellations are studied in 

^/^ . details. 

<^ . Keywords: Gibbs point process, random tessellations, stochastic geometry, pseudo-likelihood 

jrt I estimator, spatial statistics. 



1. Introduction 

In the domains of physics and biology, some large-scale random geometric structures can be 

mathematically modeled using Poisson-Voronoi or Poisson-Delaunay tessellations. In cosmology 

^O . for instance, sinc e |20l| . modeling the large-scale galaxy distribution generally relies on Voronoi 



If^ • tessellations (see [16| and [17|). In biology, Voronoi tessellations are often used to model the 

ly-T , cellular configuration of a tissue (since the seminal work of [15|). This tool is also relevant to 

(•^ ■ model the geometrical structure of proteins (cf. [2^ for a state of the art) or microstructures like 

^^ , foams. Mathematical properties of Poisson-Voronoi and Poisson-Delaunay tessellations have been 

widely studied (see [19] for instance). 

Unfortunately, these models have the disadvantage of yielding strong independence properties 
.J ■ due to the Poissonian nature of the underlying point process. In different biological studies, the 

rS I necessity to introduce an interaction between the cells of the tessellation to achieve greater realism 

j^ ' has indeed been emphasized. In [13J for instance, the interaction between neighbouring epithelial 

cells is dealt with using a Hamiltonian energy. This Hamiltonian is a function of the area of each 
cell of the Voronoi tessellation, but it involves also a pair-interaction that depends on the length 
of the common ed ge o f two cells. The same kind of interaction (but between two types of cells) is 



also considered in [12]. Moreover, some geometric hardcore interactions are sometimes demanded. 
As an example, Lautensack and Sych ([18|) modeled foams by a tessellation built from a Matern 
model with hardcore interaction. The resulting tessellation is then constrained to reach a desired 
regularity. The study of the regularity of the tessellation is also at the heart of the article of Eglen 
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and Willshaw ([11|): their work shows the relevance of forcing the geometry of cells in order to 
model retinal neurons. 

It is thus natural to consider Gibbsian modifications of the Poisson-Voronoi or Poisson-Delaunay 
tessellation, involving a smooth interaction but also a hardcore interaction (in a general sense, see 
Definition [1]), in order to produce more realistic models of interacting random structures. 

A first mathematical model has been proposed by Ord (see the discussion in [27]). In this 
model, the interaction relies on each cell of the Voronoi tessellation. In particular, a classical 
hardcore interaction forces the cells not to be too small. This model can be viewed as a nearest 
neighbour Gibbs point process and was studied in [2|. Its existence on the infinite support W^ 
is implied by the results in [3| and [4]. A Birth-Death simulation algorithm for simulating such 
nearest neighbour Gibbs processes is presented in [5|. However, tessellations involving geometric 
hardcore interactions do not generally belong to a classical theoretical framework as the previous 
one. They are in general not hereditary in the sense that, when removing a point from an allowed 
tessellation, the resulting tessellation may become forbidden (see Sectional or [10], about this 
property) . Consider for instance a generalization of the Ord process where the cells are forced not 
to be too large: this natural model is not hereditary. The existence of a Gibbs Delaunay- Voronoi 
tessellation on IR'' associated to a large class of possible non-hereditary interactions has been proved 
recently in [8| and [9|. For these processes, no simulation algorithm has been presented so far. 

From a statistical point of view, the issue is the estimation of the interaction. Assuming a 
parametric form, this may be achieved through the maximum likelihood or the pseudo-likelihood 
procedure. The maximum likelihood estimator suffers from a lack of theoretical justifications, 
except for restrictive examples of interacting point processes (see [2^ for a review) , which do not 
concern tessellation models. On the other hand, some theoretical results are available for the 
pseudo-likelihood estimator. Consistency and asymptotic normality are proved in [7| in a general 
framework including some Gibbs tessellation models, but without any hardcore interaction. A 
generalization to interactions involving a possible non-hereditary hardcore part is considered in 
[10| . In this article, the consistency of the estimation of both the hardcore part and the smooth 
part of the interaction is proved, in a setting concerning a large class of tessellation models. 

In the present article, we rewrite these theoretical results, sometimes established in an abstract 
setting, to the framework of Gibbs Delaunay- Voronoi tessellations. Moreover, some theoretical 
complements are given. In particular, a Birth-Death-Move algorithm is presented to simulate 
Gibbs tessellations with non-hereditary interactions, and a convergence result is proved. We also 
extend the recent concept of residuals introduced in [1| to the non-hereditary setting. Nevertheless, 
the aim of this article is mainly to clarify how to handle (non-hereditary) Gibbs Delaunay- Voronoi 
tessellations in practice: what kinds of models are available? How should we simulate these 
tessellations? How should we fit them to a data set and validate the fitted model? 

In the first part, the formal definition of Gibbs Delaunay- Voronoi tessellations is given. We 
restrict ourselves to tessellations on the plane R^ for simplicity. Three example models are then con- 
sidered: a non stationary crystallized triangulation model, a stationary interacting Delaunay model 
and a Voronoi tessellation model. We think that these models could be relevant for the biological 
applications cited before. Moreover, they are used throughout the article to illustrate the proposed 
methods. In Section [3J we explain how to simulate (non- hereditary) Gibbs Delaunay- Voronoi tes- 
sellations thanks to a Birth-Death-Move Metropolis-Hastings algorithm. Some simulations of the 
three above examples are presented. In Section IH we consider the estimation issue. As explained 
there, the pseudo-likelihood approach is preferred to the maximum likelihood procedure for practi- 
cal reasons. As a matter of fact, the maximum likelihood estimator is prohibitively time-consuming 
in our setting. However, if possible, maximum likelihood could be used in a second step to refine the 
pseudo-likelihood estimation. A procedure is presented to estimate both the hardcore parameters 
and the interaction, as considered in [10| . Finally, the concept of residuals as recently introduced 
in [l| is generalized, which gives a method to validate the fitted model. In the appendix, we present 
some theoretical justifications. They concern the existence of Delaunay- Voronoi tessellations, the 
convergence of the simulation algorithm, and the consistency of the estimation procedure. 



2. The Gibbs Delaunay-Voronoi tessellations model. 

2.1. The Poisson Delaunay-Voronoi tessellations. 

In paragraph l2.1.1l we recall the basic definition of point configurations. In 12.1.21 some reg- 
ularity assumptions are given to ensure that the Delaunay-Voronoi tessellations are well-defined. 
Randomness is introduced in paragraph 12.1.31 via Poisson point processes, to define the well- 
known Poisson Delaunay-Voronoi tessellations which are models of random tessellations without 
interaction between the cells. The interaction is introduced in Section 12.21 



2.1.1. Point configurations. 

Let us denote by R^ the 2-dimensional Euclidean real space. B{R^) is the set of bounded Borel 
sets in R^. The state space A^(R^) is the set of regular locally finite point configurations 7 in R^ 
defined by 

-a) for aU A in B{M?), Card(7 n A) < +00 \ 

A^(R^) = •^ 7 C R^ such that -b) four points of 7 are not on a same circle > , (1) 

-c) for every half plane iJ in R^, Card (7 n H) > J 

where Card(7 H A) denotes the number of points from 7 in the set A. Let 7 be in A^(R^) and A a 
Borel set in R^ , we denote by 7a the restriction of 7 on A which is just the set 7 n A. For a point 
X in M^, we denote by 7 4- x the configuration 7 U {a;} and if x belongs to 7, 7 — a; denotes the set 
7\{a;}. 

2.1.2. Delaunay-Voronoi tessellations. 

Let us recall the definition of Delaunay-Voronoi tessellations, which are given for example in 
[19| page 15. For a point configuration 7 in A^(M^), a set of three points T = {a;, y, z} belonging 
to 7 is a Delaunay triangle in 7 if the open circumscribed ball B{T) of T does not contain any 
point of 7. The Delaunay tessellation Del{j) is defined by the set of all Delaunay triangles T in 
7. By points b) and c) in ([1]), 'Del{j) is a partition of the space M^. 

Concerning the Voronoi tessellation coming from 7, for every a; in 7, we define the Voronoi cell 
C{x, 7) by 

C(a;,7) = |z e R^ such that Vy G 7\{a:} |z - a;| < |z - y||. 

From points a) and c) in ^, we remark that C(a;,7) is a bounded closed convex set in R^. The 
Voronoi tessellation Vor{'y) is defined by the set of all C{x, 7) for x in 7. Vor(7) is also a partition 
of the plane R^ . 

There are some relations between these two tessellations. Indeed, T = {x,y,z} in 7 is a 
Delaunay triangle if and only if C{x, 7) n C{y, 7) n C(z, 7) j^ 0. 

2.1.3. Poisson Delaunay-Voronoi tessellations. 

In this paragraph we define the Poisson Delaunay-Voronoi tessellations as in [l9| or [29|. Let 
us recall that the space of point configurations A^(R^) is endowed with the cr-algebra (t(A^(R^)) 
generated by the sets {7 e A^(R^), iVA(7) = n}, n e N, A e i3(R^), where A^a(7) denotes the 
number of points of 7 in A. The most prominent probability measures on A^(K^) are the Poisson 
processes. Let us denote them by tt", where v \s & locally finite measure on R^ and stands for the 
intensity measure (see |1^ page 83). When v is equal to z\ (z > 0, A the Lebesgue measure) we 
simply write tt^ which represents the classical stationary Poisson Point Process with intensity z. 
Let us remark that tt'^ is not necessary stationary but obviously tt^ is. 

For every A in S(R^), tt^ (respectively tt^) denotes the Poisson process tt" (respectively tt^) 
restricted on A. 

The law of Vel{j) (respectively Vor{j)) under the process tt" is called the Poisson Delaunay 
(respectively Voronoi) tessellation with intensity ly. These Poisson Delaunay-Voronoi tessellations 
are well-studied in [19?|, Section 4. 



2.2. Random Delaunay-Voronoi tessellations with interaction. 

This section is devoted to the presentation of interacting random Delaunay-Voronoi tesseha- 
tions. The interaction is introduced, via Gibbs modifications of the Poisson Delaunay-Voronoi 
tessellations, by specifying the conditional densities. This is the classical strategy used in physics 
and biology (see for example |28|). 

For every A in B{M?), we consider the conditional density /a with respect to the Poisson process 
TT^ defined by 

/A(7A,7A0 = ^7^e-^-(-^-'^-), (2) 

where 7a is a point configuration inside A, 7ac is a point configuration outside A and -Ea(7a,7a<=) 
is the energy of 7a given the outside configuration 7ac. i?A is a functional from Vor{'^) or Vellj) 
to R U {+00} which we will precise later. Z\{-/\a) :— J e~^'^^^'^'^'^''-'7rJ(((i7A) is the normalization 
constant in order to have a probability density under tt^. 

Let us remark that the conditional densities favor the point configurations with low energy and 
conversely penalize the point configurations with high energy. If the energy EAijAjjA") is equal 
to infinity then, with probability one, the configuration is even forbidden and does not appear. 
One classical example of such a situation is when the points of 7 are prohibited from being closer 
than a distance R apart, id est -Ea(7a,7A'^) = +00 if there exist x in 7a and y in 7 such that 
|x — 2/1 < R. This constraint is usually designated as a hardcore interaction. In this paper, we 
generally call hardcore interaction any situation where -Ea(7a,7a<=) = +00. 

Definition 1. An energy (or an interaction) is said to contain a hardcore part if EA^'yAj^A") = 
+00 for some 7 and some A. 

We denote by A^oo(R^) the set of allowed configurations which is defined by 

Xoc(R^) = {7eX(K2) suchthat for all A inS(R2), ^^(^^^^^,) < +oo|. (3) 

Now let us define the model of random Delaunay-Voronoi tessellations with interaction. 

Definition 2. A probability measure P on A^oo(R^) is a Gibbs Delaunay-Voronoi tessellation for 
the energies {Ea)agB(R^) '^'^'^ ^'^s intensity measure v if for every A in BiJS?) and for P-almost 
every outside configuration ^A", the law of P given ^A" is absolutely continuous with respect to tt"^ 
with the density /a(-,7a<^)- 

This definition of Gibbs measures is the classical one that can be found for example in |25| . 
Let us point out several problems about the existence of these Gibbs Delaunay-Voronoi tessel- 
lations. First of all there are some conditions on the energies (-EA)Ae6(K2) to ensure that the 
conditional densities (/a)agB(k^) ^-re well-defined and compatible. Moreover, even if it is the case, 
it is not obvious that Gibbs Delaunay-Voronoi tessellations exist and are unique (the non unicity 
of Gibbs processes is called phase transition in statistical mechanics). In Section 5.1, we give some 
conditions to ensure the existence of stationary Gibbs Delaunay-Voronoi tessellations for energy 
functions given by (j4]) and ([5j below. As far as we know, uniqueness or phase transition for such 
Gibbs measures have not been proved. For similar Gibbs models with multi-type particles, a phase 
transition result is given in |6|. 

Since we are interested by models of random Delaunay-Voronoi tessellations, the energy func- 
tions have to depend on the local geometry of the Delaunay or Voronoi tessellations. We need 
some notations. Let us first define the neighbour relations ^-Dei, ^Vor between the cells in 'Del{'y) 
or Vor ('-/). 

For aU T, T' e VeKcf), T ^vei T' if Card(T n T') > 2. 
For aU C, C E Vor{-f), C ~vor C" if C D C ^ 0. 

In fact, T ^-Dei T' if T and T' have a common edge and C ^vor C if C and C" have a common 
edge at their boundary. 



Now let us define the cells in Vell^) or Vor('-/) which are inside or outside a given bounded 
set A in S(R^). A triangle T e 'Del{'-)) (respectively a cell C G Vor(7)) is outside A if for every 
configuration 7^ in A, T (respectively C) is in 'Del{j\': U 7^) (respectively Vor(7A<: U 7^))- In- 
other words, T (or C) is outside A if T (or C) remains in Vell^) (or Vor{^)) for any modifica- 
tion of the configuration 7 inside A. T (or C) is inside A if it is not outside A. We denote by 
'DelAi'y) (respectively VorAi'y)) the cells T in 'Del{'-f) (respectively C in Vor{'-f)) which are inside A. 

a) A general form for the energy of a Delaunay tessellation 
We define the energy Ea of the Delaunay tessellation by 

i?A(7A,7A0= "E ^i(r)+ y. V2iT,r), (4) 
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T' in DcIa 


(7) 



where Vi is a function from the space of triangles T to M U {+cx)} and V2 is a symmetric function 
from T^ to M U {+00}. In Section 2.3, we give precise examples of functions Vi and ¥2- 

b) A general form for the energy of a Voronoi tessellation 
Similarly, we define the energy Ea of the Voronoi tessellation by 

i^A(7A,7A0= E ^i(^)+ E V2{C,C'), (5) 

CGVorA(7) {C,C'}CVor(7) 

C or C' in VorA(7) 

where Vi is a function from the space of bounded convex sets C to MU {+00} and V2 is a symmetric 
function from C^ to MU{+oo}. In Section 2.3, a precise example of functions Vi and V2 is provided. 

2.3. Three explicit reference models 

We present in this section three explicit examples of Gibbs Delaunay- Voronoi models, that will 
be our reference models until the end of the paper. All the functions Vi and V2 given in this section 
satisfy the assumptions in Section 5.1 and so the associated Gibbs Delaunay- Voronoi tessellations 
exist. 

Model 1: a non stationary crystallized triangulations model. 
In this model, we propose to define a non stationary random triangulation in which the angles of 
triangles are forced to be larger than a fixed real a in [0, ^[. If a is chosen close to ^, the model 
produces rigid random triangulations. Moreover, it is possible to have a non stationary density of 
points. 

We assume that the intensity measure v is absolutely continuous with respect to the Lebesgue 
measure A and the energy Ea is defined by ^ with 

^i(^)=| otherwise, ^'^^ ^^ = 0, 

where a{T) is the minimal angle inside T . In fact, the energy Ea is equal to plus infinity if there 
exists a triangle inside A which is too fiat. 

Model 2: a stationary interacting Delaunay model. 

In this model, we propose to study an example of stationary Delaunay triangulation with 
interaction. It is a simple model in order to present different practical and theoretical aspects 
in this work (modeling, simulation, estimation). This model has not the ambition to be directly 
usable in physics or biology. 



First of all. we fix the intensity measure v to be equal to z\. Via a geometric hardcore 
interaction, we force the edges not to be too small, the triangles not to be too large and via a 
smooth interaction the large perimeters of triangles are favored or penalized (depending on the 
sign of Q). More precisely, let < e < a and Q be in M. the energy Ep^ is defined by ^ with 

r +00 \il{T)<e, 

Vx{T)=\ +00 if i?(T)> a, and ^2=0, 

[ dVer{T) otherwise, 

where 1{T) is the minimal length of sides of T, R{T) is the radius of the circumscribed ball of T 
and Ver{T) is the perimeter of T. 

Model 3: a Voronoi tessellation model. 

In this third example, the interaction is defined as far as possible to fit with the biological 
applications evoked in the introduction (geometric regularities of cells, interaction between cells), 
although other interactions could be chosen. We suppose that the model is stationary so we fix 
the intensity 1/ to be equal to zX. The geometry of cells is controlled by a hardcore interaction 
Vi which forces the cells not to be too small, too large or too flat. Moreover, we add a smooth 
interaction involving a competition between the volumes of neighbours cells. 
Let < £ < a, S > 1/(2V^) and 9 be in R. The energy Ea is defined by ^ with 



Vi ■.C^Vi{C) = 



FOO if /imin(C) < e, 

FOO if /imax(C) > a, 

Fcx) ifhl^,iC)>BVol{C), 

otherwise. 



V. : (CO ^ V.iCC) ^ 9 ('":^l!i;^.\^1£l^ - 1 



and 

TA • (n r'\ ^ v^(r r"\ = a I — 

j'nm{Vol{C),Vol{C')) 

where /iinin(C') is the minimal distance between the center x of the cell C and the edges of the 
boundary of C. Similarly, /imax(C') is the maximal distance between x and the edges of C (see 
Figure [J) . 




Figure 1: Example of Voronoi cell C with center x and distances h^ 

The choice of the power i is arbitrary and may be changed. Nevertheless, it seems to lead to 
more realistic simulations. The parameter B controls the form of the cell: the smaller B, the more 
regular the cell. For instance, for a hexagonal cell, B — l/(2-v/3) ~ 0.29. Let us remark that if 9 
is positive the interaction V2 forces the neighbour cells to have the same volume. Conversely, if 9 
is negative it forces the neighbour cells to have different volumes. The sign of 9 is crucial in this 
model. 



3. Simulations 

3.1. Gibbs tessellations on a finite window 

In this section, we deal with the simulation of our models. First of all, let us remark that Gibbs 
Delaunay-Voronoi tessellations are processes on R^ and so one has to restrict or approximate them 
inside a fixed window A. The most natural choice would be to simulate the restriction on A but 
it is almost impossible to do it. Therefore the common method is to consider the finite volume 
Gibbs approximations on A, which is the probability measure absolutely continuous with respect 
to the Poisson process 7rJ( with the density /a given in ^. The outside point configuration 7a <= 
is then fixed and chosen arbitrarily. Results in statistical mechanics show, in general, that the 
thermodynamic limits of these finite volume Gibbs measures, when the volume A goes to R^, are 
Gibbs measures (see [25|). Therefore, finite volume Gibbs measures are good approximations of 
our models. 

There are essentially three possibilities to fix the outside configurations. The first possibility is 
to consider the empty outside configuration but it is not usable in our context because it produces 
non bounded Delaunay-Voronoi cells and so non computable energies. The second possibility is to 
fix an outside point configuration which has to be specified explicitly. Again, this is not practical 
because the strong hardcore interaction, which appears in our three reference models, makes it 
difficult to find such a configuration. The third possibility is the periodic outside configuration 
which is built by periodization in the full plane R^ of the random configuration inside A. We 
choose this approach because it seems relevant to deal with the hardcore problems coming from 
the boundary effects. 

Now, let us give the precise construction of the periodic finite volume Gibbs measure. The 
simulation window A is chosen as the square [0, 1]^. Any other window in R^ can be considered 
in the same way. A simple rescaling procedure enables us to reduce to this case. For every point 
configuration 7 in [0, 1]^, we denote by 7 the periodic configuration on R^ defined by 

where Ti is the translation in R^ with respect to the vector i. 

In the sums over the cells in (jl]) or ([5]) , we must ensure that each collection of periodic cells has 
a unique contribution in the computation of the periodic energies. A solution consists in selecting 
only the cells whom barycenters are in [0, 1]^. Therefore, for any Voronoi cells C (respectively 
Delaunay triangle T), we denote by < C > (respectively < T >) the barycenter of the cell C 
(respectively triangle T). Similarly, for any couple of Voronoi cells {C,C') (respectively trian- 
gles {T,T')) we denote by < C,C' > (respectively < T,T' >) the barycenter of the set C U C" 
(respectively TUT') . 

The periodic energy E{j) associated to the energy of a Delaunay tessellation (j4]) is defined by 

^(7)= E ^i(^)+ E V2iT,r). (7) 



TeVelij) 


{T,T'}CVel{j) 


<T>e[0,l]^ 


C~-DslC 




<T,T'>e[0,l]^ 



Similarly the periodic energy E{'~f) associated to the energy of a Voronoi tessellation ^ is 
defined by 

^(7)= E ^i(^)+ E y2ic,c'). (8) 

CeVor(7) {C,C'}CVor(7) 

<C>e[0,l]^ C~vorC' 

<c,c'>e[o,i]^ 
Now let us define the Gibbs process on [0, 1]'^ with periodic outside configuration. 



Definition 3. The periodic Gibbs Delaunay-Voronoi tessellation P is the point process on [0, 1]^ 
which is absolutely continuous with respect to the Poisson point process on [0, 1]^ (denoted by tTq), 
where the density f is defined, for every 7 in [0, 1]^, by 

f{-f) = ^e-^'--'^ and Z = f e-^'^^'K^{d-/'). 

Remark 1. In the case where the intensity measure and the energy functions are stationary, we 
know there exists at least one stationary Gibbs measure (see Theorems 1 and 2 in the appendix). 
However, non stationary Gibbs measures may exist too. This phenomenon is called the breakdown 
of symmetry (see ]2a] 4-1, for instance). It is not proved theoretically that symmetry breakdown 
can occur for the models studied in this paper. Sometimes, this phenomenon can be observed via 
simulations. Nevertheless let us remark that this will not be possible with the periodic Gibbs models 
simulated here. Indeed, if the intensity measure and the energy functions are stationary, as in 
models 2 and 3, then the associated periodic Gibbs models are stationary too and the symmetry 
breakdown is not observable. Free or configurational boundary conditions should be investigated. 

3.2. Algorithm of simulation 

As explained above, to deal with the boundary effects, we choose to simulate periodic tessella- 
tions. For sake of brevity, we confuse 7 and 7 in this section, similarly we will use / instead of /. 
The window of simulation is A = [0, 1]^ and we omit the indexation by A in the sequel. Moreover, 
the notation A^oo(K^) is abusively extended in this section to the periodic tessellations with finite 
energies. 

There exist different algorithms to simulate finite volume Gibbs processes. Some perfect simu- 
lation algorithms have been developed (see [14], [26]) but they seem not to be really implementable 
in our context due to the strong rigidity of our hardcore models. So we make the choice to simulate 
them via the classical Birth-Death- Move Metropolis-Hastings algorithm, which we recall below (see 
also ^). 

For X e [0, 1]^, N{x, cr"^) denotes the Gaussian distribution on R^ centered at x with covariance 
matrix diag{a^ , cr^), where cr > 0. This law is the proposal density for moving a point. Note that 
if the moved point falls outside [0, l]'^ (in step 5. below), it is replaced inside [0, 1]^ by the periodic 
property. We assume that the intensity measure v is absolutely continuous with respect to the 
Lebesgue measure, and we denote by g its density. In particular, in the stationary case, i.e. when 
V = z\ with z > 0, then g is identically equal to z. 

1. For 70 e A^oo(R^), let n = card{^o). 

2. Draw independently a and b uniformly on [0, 1]. 

3. If a < 1/3 then generate x uniformly on [0, 1]^ and set 

71 = P° ^ ("+l)/(7o) ' 

1 7o otherwise. 

4. If a > 2/3 then generate x uniformly on 70 and set 

1 7o otherwise. 

5. If 1/3 < a < 2/3 then generate x uniformly on 70, generate y ^ Af{x, a'^) and set 

\jo-x + y if 6</(2-£+J£), 



^1 1 , • 

1 70 otherwise. 



/(70) 



6. Iterate from 1. where 70 -^ 7i- 



This algorithm can be refined: The probability of move, birth or death proposals may differ 
and may depend on x. similarly the law for choosing a point before killing it, adding it or moving it 
may be chosen properly (e.g. according to the intensity law v). The idea of this procedure is that, 
starting from an allowed configuration 70, the iterations converge to the realization of an invariant 
measure which is the Gibbs process we want to simulate. For classical Gibbs point processes, this 
convergence is proved for example in Section 7.3 of [22]. In our setting, the convergence is not obvi- 
ous. It mainly relies on the following connectivity property (see also Definition [5] in the appendix): 
from any allowed configuration 7, it is possible to reach another allowed configuration 7' thanks to 
an iterative birth-death-move procedure as above. Since a hardcore Gibbs tessellation may be very 
rigid, this property does not always hold (consider for instance the Delaunay tessellations where 
all the triangles are imposed to be almost equilateral). Yet, in most situations, the connectivity 
exists. Let us note that the moving step is crucial here, because it allows the connectivity of rigid 
tessellations that a simple birth-death procedure would not. 

We show in the appendix that, under some assumptions, the algorithm converges. These 
assumptions are fulfilled for Model 2 when 2e < a < i. For the Voronoi tessellation presented 
in Model 3, we claim that the convergence holds for a large reasonable set of parameters (e, a, 6). 
However, in this case theoretical justifications become tedious and are not achieved in this paper. 

3.S. Practical implementation 

In the above algorithm, the choice of the initial configuration is crucial. We must start from 
an allowed configuration 70 in [0,1]^, i.e. 70 € A^oo(K^). For this reason, we cannot start from 
the empty configuration. In our simulations, we chose to start from the point configuration whom 
Delaunay tessellation is a regular lattice of triangles (see its plot on top left of Figure [5]) . This 
starting configuration is allowed by all our hardcore models, provided the distance between points 
is properly chosen. 

The computation of the ratio in steps 3.-5. of the algorithm is time-consuming since it supposes 
the computation of two tessellations plus a calculus on their cells to obtain their respective density. 
But it is possible to simplify this computation by focusing on a smaller window. Indeed, consider 
for instance step 3, i.e. the birth case (the same approach remains true for the death or the move 
step). When one adds a point x in a configuration, the new tessellation differs from the previous 
one only in a neighbourhood of x. Thus, the ratio of the densities reduces to a difference of energies 
in this neighbourhood. The size of this neighbourhood is determined by the size of cells around 
X. For instance in Model 2, the diameter of cells is forced to be smaller than 2a, so it suffices to 
focus on a ball with radius 2a around x to compute the difference of energies. 

When there is a hardcore interaction, the convergence of the algorithm may be slow. Indeed, 
when one adds a point, the new tessellation can be forbidden. Hence, in presence of a strong 
hardcore interaction, most of the new tessellations proposed by the birth or death step in the 
algorithm will be refused. For this reason, we check the progress of the algorithm by monitoring 
control. Every one thousand iterations, we count the total number of points in the configuration 
as well as the number of accepted birth steps (among one thousand steps), the number of accepted 
death steps and the number of accepted move steps. The plot of theses numbers all along the 
iterative algorithm helps us to check the stabilization of the iterative process (though this is not a 
proof of the convergence to the invariant measure) . 

3.4- Some examples 

We present some simulations of the models introduced in Section [^31 We use the Birth-Death- 
Move algorithm presented before with a = 0.015. 

Figure [2] shows a simulated tessellation from Model 1 where a — it/&. Let us denote by g the 
density of the non-stationary intensity measure v. It is 1-periodic with respect to each component, 
i.e. for aU {x,y) G R^^ ^(3,^ ^,2/) = g{x,y) and g{x,y + 1) = g{x,y), and for {x,y) € [0, l]^ we 
assume 

5(x, y) = 100[(x - 0.5)2 + (^ _ q, 5)2]-0-75. 



The triangles of such a tesseUation are forced not to be too flat and to be more dense around the 
point (0.5,0.5). The monitoring control seems to justify the convergence of the algorithm after 
5. 10** iterations. 

In Figure 131 two tessellations from Model 2 have been simulated with a = 0.08, z = 1000, and 
9 = ±5. We did not introduce in these simulations the hardcore parameter e. When 9 > (on 
bottom), the tessellation is more likely to exhibit a small number of vertices since the total sum 
of perimeters is then low, which minimizes the energy. For 9 < (on top), this is the contrary: 
the energy is minimal when the total sum of perimeters is high, inducing a lot of vertices in the 
tessellation. In these two cases, the size of the triangles is controlled by the hardcore parameter 
a, which might be unnecessary when 9 < but is certainly a big constraint when ^ > 0. 

Model 3, involving Voronoi tessellation, has been simulated for a = 0.05, B = 0.625, z = 100, 
and 9 = —0.8, —0.5, 0.5, 0.8. The hardcore parameters a and B force the cell not to be too large 
or too flat. We did not impose a minimal length of sides through the hardcore parameter e. The 
parameter 9 quantifies the dependence between the size of neighbour cells: when 9 < they are 
most likely to have different sizes, whereas for 9 > they tend to exhibit the same volume. These 
two opposite behaviors are clearly observable in the two extreme cases 9 ~ —0.8 and 9 — 0.8 in 
Figure m When \9\ = 0.5, this difference is more difficult to distinguish. A challenging task will 
be to properly estimate the parameters a, B, z, and 9 in both the apparently closed situation 
6 = —0.5 and 9 = 0.5. This problem is addressed in the next section. 

In Figure [SJ one can check the convergence of the algorithm for the four above simulations. 
Note that in the very rigid case 9 — 0.8, the iteration process needs a lot of time (tens of thousands 
iterations) before starting a birth-death step. This shows the importance of the moving step. 

More simulations of Model 3 are presented in Figure [HI with their monitoring control in Figure 
[71 They show the impact of B on the geometry of the cells. It plays a bigger role when 9 < 0, 
since in this case we note the cells may be very flat without this hardcore parameter {B = +oo). 
Let us remark that some clusters of small cells appear when 9 < 0. This is more visible when B is 
not too constraining. 
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Figure 2; Simulation of Model 1 with a = tv/6 and g{x,y) = 100[(a;-0.5)^ + (y-0.5)^]"'''^'"' (5.10^ iterations). 
Top left: initial point configuration; Top right: monitoring control (from top to bottom: number of moved 
points, number of birth points, number of killed points and total number of points, pointed out every 1000 
iterations); Bottom: final simulated tessellation. 



11 




Nb of points 




■■\1\/W\.\fv vA-vV^-V>. 




Iterations ('1000) 



^^l^A 


^ 


(lUMvi/\l 


'^HaJ^w^^^ 





50 


100 

Births 


150 200 


'jl'Af^Aj' 


^^vv|vA 


•iHfA^fvWV 


^^fMlVA^V 





50 


100 

Kills 


150 200 


M|vylAVVv^A4^\|^^^ 





50 


100 

Nb of points 


150 200 


AyH_V..,..Av/-v-MA-^^ 



Figure 3: Simulation of Model 2 with a = 0.08, z = 1000 and 9 = -5 (top), e = 5 (bottom) after 2.10^ 
iterations. Top left: final simulated tessellation when 8 — —5; Top right: monitoring control when 9 = —5 
(with the same plots as for Figure[2|; Bottom left: final simulated tessellation when 9 = 5; Bottom right: 
monitoring control when 9 = 5 (with the same plots as before). 
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Figure 4: Simulation of Model 3 with a = 0.05, B = 0.625, z = 100 and 6 = -0.8 (top left), 6 = -0.5 (top 
right), 6 = 0.5 (bottom left) and 6 = 0.8 (bottom right). These are the final simulated tessellations after 
2.10^ iterations when \6\ = 0.5 and 5.10* iterations when |^| — 0.8 (see the monitoring control in Figure 



13 



^W■'%^,^V-^^vA.^tw^^ ^ ^•^''^/\k||^^^^|]|AfY^^ 



100 200 300 400 500 



'^VWs 



'vl-/WJ^/V--'Wv/^VvvvV'-'^A"v/-*v^'^^ 




100 200 300 400 



100 200 300 400 500 



■--^^xrKfA/'r^-^''^.r^K^^^A^ 



Hl\ij\f|^4/^fy^^ ^- PV^/Vvii*^H^ 



Nb of points 




100 200 300 400 

Iterations (*1 000) 



Nb of points 




Iterations ('1 000) 



iVvl\A^'^ivfViVi|uf^^ 



100 200 300 400 500 



yikikwk^^^^^^^^^^^^^ 



100 200 300 400 500 



'Li.iAl^l/LmmLAffiLA^AVjiA/r''VUA^^ 



100 200 300 400 500 



Nb of points 




Iterations (*1 000) 



Nb of points 



100 200 300 400 

Iterations (*1 000) 



Figure 5: Monitoring control for the simulations of Model 3 presented in Figure [H in the same order (from 
top left to bottom right: 9 = —0.8, —0.5, 0.5, 0.8). They consist in the same plots as for Figure (2] 
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Figure 6: Simulation of Model 3 with a = 0.05, z = 100, B = +(x (left), B = 1 (right), 9 = -0.5 (top), 
9 = 0.5 (bottom). These are the final simulated tessellations after 1.5.10'' iterations (see the monitoring 
control in Figure [7|). 
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Figure 7: Monitoring control for the simulations of Model 3 presented in Figure [6l in the same order 
{B — +00 (left), 5 = 1 (right), d = —0.5 (top), 6 = 0.5 (bottom)). They consist in the same plots as for 
Figure [1 



16 



4. Estimation 

In this section, we focus on stationary Gibbs tessellations, in particular we suppose that the 
underlying Poisson process is tt^, the stationary Poisson point Process with intensity zX where 
A is the Lebesgue measure. This is the case of Model 2 and 3 presented before. We apply our 
estimation procedure to the parameters z, e, a, B and involved in these models. An interesting 
generalization would be to include the estimation of a non-stationary intensity (as in Model 1) to 
the estimation of the parameters of the interaction. We do not deal with this task in this paper. 

Let us specify some notations. We assume that we observe a tessellation coming from a point 
configuration 7 on a window A„ = [— n, n]^. This tessellation is defined through an energy function 
£'a(7A)7a<=) as in ((3]) or (O. In the following, we denote it by S^' (7a,7a<:), since we assume a 
parametric form. We actually suppose that it depends on parameters (/3,6'): /? is the hardcore 
parameter (parameterizing the finiteness of the energy function), while 6 is the smooth interaction 
parameter. For instance, in Model 2: P — (e,a), 6* = 6*; in Model 3: /? = (e,a,i3), 9 = 9. This 
distinction is presented clearly in |AppendixA.3[ where some theoretical results for the asymptotic 
consistency of our estimators are given. 

We consider a classical two-step estimation procedure. We first estimate /3, then we estimate 9 
and z by pseudo-likelihood where /3 is replaced by its estimator. 

The choice of the pseudo-likelihood approach instead of the classical maximum likelihood esti- 
mator is mainly imposed by practical reasons. Indeed, maximum likelihood requires the estimation, 
by simulations, of an unknown normalizing constant. This approach demands to simulate several 
tessellations according to the model, which is extremely time-consuming in the situation when a 
hardcore interaction is involved (see previous section). Moreover, the pseudo- likelihood procedure 
has the advantage of being asymptotically consistent for a large class of models (see [lO|), which has 
not been proved for the maximum likelihood estimator in such a general setting. However, when 
the hardcore interaction is not too strong, the maximum likelihood estimation may constitute a 
second step to refine a pseudo-likelihood approach. 

There is a major difficulty to overcome in order to implement the pseudo-likelihood estimation 
in our case: the hardcore interactions are not necessarily hereditary. An interaction is hereditary 
if, for every forbidden point pattern 7, then, for every point x, the configuration 7 + a; remains 
forbidden. This is equivalent to: for every allowed point configuration 7, then for every point 
X S 7, the configuration 7 — 2; remains allowed. In other words, an interaction is hereditary 
if one can remove any point from 7. This property concerns only the hardcore interaction. So 
every interaction involving no hardcore part is necessarily hereditary. The models presented in 
Section 12.31 are not hereditary. Indeed, if one removes a point from an allowed tessellation, the 
new tessellation may contain cells that are too large (for instance). As a consequence, we must 
modify the classical pseudo-likelihood contrast to take into account the so-called removable points 
as introduced in LlDi] (see Definition \^ . 

4.I. The two-step procedure 

The first step consists in estimating the hardcore parameter /3. Let us first assume that /? is a 
one-dimensional parameter. We suppose the following inclusion 

if /? < /3' then VA, E^/{-,f,,-/Ac) < +^ ^ ^f ''(7a,7aO < +^- (9) 

In this case, a consistent estimator of /3 is 

/3 = inf{/3 > 0, £;^f (7A„,7AfJ < +^}. (10) 

If instead of (j9]), the converse implication holds, then it suffices to replace the infimum by a 
supremum in (fTO|). 

In the case of a multi-dimensional hardcore parameter /3, we estimate each of its components 
as above. 



17 



For instance, for Models 2 and 3 presented in Section 12. 3[ Property ([9]) is satisfied by the 
hardcore parameters a and i?, while the converse holds for e. As a consequence, following (jlOl) . 
natural estimators for these examples are (the notations are the same as in[ 



For Model 2: 



e = min{l{T), T e VdA^{-f)}, 
a — max{R{T), T £ I?e/A„(7)}. 



For Model 3: 



e = min{h^iniC), C G VorA„(7)}, 
a = max{/imax(C), C eVorA^ij)}, 
B = max{hl^^^{C)/Vol{C), C £ VorA^il)}- 

The second step consists in estimating the smooth interaction parameter 9 and the intensity 
parameter z. We use the pseudo-likelihood procedure for the reasons explained before. To deal 
with the non-hereditary problem, we must introduce the concept of removable points. 

Definition 4. Let 7 be in A^(R^) and x be a point of -f, then x is removable from 7 if there exists 
A e S(M^) such that x e A and 

El'^ijA- X,JA'=) < +00. (11) 

The following proposition, proved in [lOj . gives a more intuitive approach and justifies the name 
of removable points. 

Proposition 1. Let 7 be in A^oo(K^) and x be a point of j, then x is removable from 7 if and 
only if J — X is in Mac (K^ ) • 

From the definition, it is clear that the property of being a removable point from 7 depends only 
on the hardcore parameter (3 and not on 6* or z. Thus, we denote by 7?.^ (7) the set of removable 
points in 7. 

The more rigid the tessellation, the less removable points there are. In particular, if there is no 
hardcore part in the interaction function, every point of 7 is removable. In Figure[51 the removable 
points of previous simulations are encircled. 

We are now in position to introduce the pseudo-likelihood contrast function, adapted to the 
non-hereditary case: 

PLLa„(7,z,/3,0) = / zexp{-h^^\x,j))dx+ ^ (/i'^'^(x,7 - x) - ln(z)), 

where h^'^{x, 7 — x) is the local energy of x in 7, defined for every x £ TZ^{'y) by 

h^'\x,j- x) = E^/{^A,1A^) - Ei-'ilA - x,7aO, (12) 

where A is a set containing x as in Definitional Let us point out that h"' {x,j) is just equal to 
h^'^{x, {j + x) — x) and is always well-defined for 7 in A^oo(K^)- 

The parameters 6 and z are estimated by minimizing PLLa„-, where the hardcore parameter 
/? is replaced by its estimator /3 obtained in the first step: 

(z, 6) = argminz^gPLLA„{j,z,(3,0). (13) 

The consistency of this estimation procedure is considered in |AppendixA.3[ 
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Figure 8: Removable points (encircled) from: The Delaunay tessellation simulated in Figure |3] where 9 = 5 
(left); The Voronoi tessellation simulated in Figure |4] where 6 = —0.5 (right). 



4-2. Practical implementation 

The opthnization oi PLL\^ requires the calculus of the local energy h^'^[x,^) for any x G A„. 
This is the same calculus as the one needed in step 3 of the algorithm presented in Section |3^ and, 
as explained there, it can be achieved by focusing on a window around x. Moreover, this compu- 
tation requires the knowledge of 7a<:, the configuration outside this window. To prevent boundary 
problems, it is actually necessary to compute PLL on a sub-window of the initial observation 
window A„. We denote abusively A„ this sub-window in the following. 

The derivative of PLL\^ with respect to z yields the following estimator for z: 



/^^exp(-/i'^.S(x,7))da: 
card{n^ (j) n An) 



(14) 



where card{TZ^ {'y) n A„) is the number of removable points from the observed point pattern 7 in 
A„. 

In the simple case where is a one-dimensional parameter and h°'^{x, 7) is a sufficiently regular 
function, the minimization of PLLa„ in 9 can be reduced to the determination of the root of an 
equation. Indeed, from the derivative of PLL\^ with respect to 9, we obtain in this case that 9 is 
the solution of 



dhP 

, ~de' 



■(x,7)exp(^-/i'^'^(a;,7) 



dx 



E 



a;eK'5(7)nA„ 



~de' 



■{x,j-x), 



(15) 



where A;^ = {a; G A„, j + x e A^oo(K^)}- 

Moreover, when h^'^{x,j) depends linearly on 9 (as in Models 2 and 3), for all x such that 

j + xeMc^{M.^), 

-^^{x,j)^9h^-Hx,j), 

which simplifies equation (J15p above. 

From a practical point of view, we first estimate 9 thanks to (|15p . where z is replaced by z 
given by (HH). Then we deduce z by plugging 9 into ^^. In both these estimations, the involved 
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integrals are approximated by Monte Carlo (this is the most time-consuming step of the estimation 
procedure). 



4-. 3. Some examples 

4.3.1. For Model 2 

We implement the estimation procedure on simulations of Model 2. We do not introduce the 
hardcore parameter e here. The estimation of a, 9 and z has been done from 200 replications of 
Model 2 when a — 0.08, z — 1000 and 9 — ±5, simulated as in Section [X^ The results are shown 
in Figure [9] and [101 We have distinguished two cases: first estimating 9 by supposing z = 1000 
known, then estimating both 9 and z. In this last case, one can note in the bottom right plot of 
these figures the closed relation between z and 9. Although the models that we consider are well 
identifiable, it is not surprising to observe this closed relation: it is implied by the Euler's formula, 
which connects linearly the number of cells and the number of vertices in a tessellation (see 3.2.11 
in [19(1). Therefore, if z is under-estimated, 9 will tend to be under-estimated as well, in order to 
respect this linear relation. 

When 9 = —5 (FigurelHl), the simulated tessellations rely on about 1500 points and all of them 
are removable. The estimation of a when 9 = —5 actually shows that this hardcore parameter 
is useless in this case: the cells of the tessellation naturally satisfy the hardcore condition. It is 
interesting to note that this misspecification does not affect the estimation of the smooth interaction 
parameter 9. The average of 9 is about —5, while its standard deviation is 0.4 when z is known 
and 1.6 when z is estimated. The average of z is 1002 and its standard deviation 145. 
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9 when z is estimated z Scatterplot of (6, z) 

Figure 9: Estimation of Model 2 when a — 0.08, 9 = —5, z = 1000, from 200 replications. 

When 9 = 5 (Figure [TO)) . the hardcore plays an important role in the model. It is well estimated 
with a standard deviation of 3.10"^. The standard deviation of 9 is 0.3 when z is known and 1.9 
when z is estimated. The average of z is 1049 and its standard deviation 313. These estimations 
seem less accurate than when 9 = —5. This certainly comes from the fact that, when 9 = 5, our 
simulated tessellations on [0, 1] x [0, 1] rely only on 500 points. Most of these points are removable 
(more than 90%), as showed in the left example of Figure [H 
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Figure 10: Estimation of Model 2 when a = 0.08, 9 ='a, z = 1000, from 200 replications. 



4.3.2. For Model 3 

Two hundred replications of Model 3 where a = 0.05, B — 0.625, z = 100 and — ±0.5 have 
been simulated according to the algorithm presented in Section [X^ (see Figure 0] for an example). 
As above, the hardcore parameter e was not introduced here. The results of the estimations are 
shown in Figure [TT] when 6 = —0.5 and in Figure [T^ when 9 = 0.5. Two situations are considered, 
assuming z — 100 is known or not. The particularity of these simulated Voronoi tessellations is 
their rigidity. The hardcore interactions are strong, forcing the cells not to be too large (through 
a) neither too flat (through B). This is confirmed by the accuracy of their estimation in both 
cases — ±0.5 (see the histograms in Figures [TT] and [T^ . But, as a consequence, there are only a 
few removable points, making the estimation of the smooth interaction parameters more difficult. 
Yet, it appears from these simulations that, in spite of the apparent similarity of the tessellations 
when 6 — —0.5 and Q — 0.5 (see Figure 2]) and in spite of the few number of removable points, the 
estimation procedure is mostly available to properly distinguish them. 

When Q = —0.5 (Figure [TT|) . there are in average 45 removable points on 265 points. The 
estimation of d remains correct: the average and the standard deviation of Q are respectively 
—0.52 and 6.4 10"^ when z is known, and —0.56 and 14.5 10"^ when z is estimated. The average 
of z is 94 while its standard deviation is 45. 

When B — 0.5 (Figure [T^ . there are only 3.5 removable points in average on about 215 points. 
In this latest extreme case, some estimations of Q and z were even impossible since there were no 
removable points at all (in 5 percent of the simulations). This shows the limit of the estimation 
procedure in presence of a very rigid tessellation. The average of d is 0.55 and its standard deviation 
is 22.8 10"^ when z is not estimated. This is surprisingly reasonable in view of the few numbers 
of removable points. When both z and d are estimated, the results become bad: the average of Q 
is 0.53 with a standard deviation of 48 10~^ and the average of z is 189 with a standard deviation 
of 345. Their joint distribution is plotted on bottom left of Figure [T^ A zoom in is plotted on 
bottom middle, where more than 90% of the points are remaining. The last plot on bottom right 
shows the repartition of Q according to the number of removable points in the tessellation, when 
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Figure 11: Estimation of Model 3 when a = 0.05, B = 0.625, 6 = -0.5, z = 100, from 200 replications. 

z = 100 is assumed to be known. There is a clear bias when the number of removable points 
is low. Since z = 100 is fixed, this low number of removable points is associated with a strong 
rigidity, so 6 is most likely to be high. Moreover, the standard deviation of 9 seems to decrease 
with the number of removable points. This is confirmed by Table [T] which contains, for a fixed 
number of removable points card{TZ^ {'y)), the number of tessellations from our simulations having 
this number of removable points (named replications) and the standard deviation of 9n calculated 
from these tessellations (denoted sd{9)). 



card{n^{-/)) 


1 


2 


3 


4 


5 


6 


7 


>7 


sd{9) {xlO-2) 


27.6 


14.9 


18.4 


15.1 


9.6 


10.3 


8.8 


8.7 


replications 


23 


41 


41 


46 


29 


19 


8 


7 



Table 1: Standard deviation of 9 according to the number of removable points, from replications of Model 
3 with 9 = 0.5. 



4-.4- Analysis of residuals 

When fitting a model to a data set, the analysis of the residuals is a standard way to check 
the quality of the model. The concept of residuals for spatial point processes is not simple. A 
general definition is proposed in [lt|, where the authors also present several diagnostic tools based 
on residuals. The definition relies on the Campbell equilibrium equation due to Nguyen and Zessin 
(see [23]), where the Papangelou conditional intensity is involved. In our context, the Papangelou 
conditional intensity does not always exist, because the hardcore interactions are not necessarily 
hereditary (see Remark 2 in [lO|). Yet, a Campbell equilibrium equation still holds, provided we 
restrict the support to the set of removable points. 
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Figure 12: Estimation of Model 3 when a = 0.05, B = 0.625, 9 = 0.5, z = 100, from 200 replications. On 
top left: A typical tessellation which is estimated. Bottom right: Repartition of 6 according to the number 
of removable points observed in the tessellation (see also Table [l]). 
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Proposition 2. Let P he a stationary Gibbs Delaunay-Voronoi tessellation as defined in Definition 
[^ For every bounded nan negative measurable function ijj from IS? x A^(R^) to M, we have 



Ef 






i^ix, ^-x)]=Ep (j^ ^{x, ^)e-''''"^-'^^v{dx)\ , 



where h^'^ is defined in il^) and Ep denotes the expectation under P . 

This proposition is proved in jlQ) . From this equation, following |l||, we can define the innovation 
process, for any bounded set A in M^ and for every function -0 as above: 






-h»'\xa) 



v(dx\ 



The residuals are then defined as an estimation of the innovations: 



R 



(a,V,/i'^'V) = Y. V'(a;,7-a^) 

xe7?.'5(7)nA 



V'(a;,7)e 



'hf>-\x,-f) 



i>{dx), 



where i> is an estimation of the intensity measure v, which is simply zX in the stationary case. 

This generalization of the residuals to the setting of possible non-hereditary interactions allows 
to perform several diagnostic plots. We refer to [l| for a presentation of different relevant choices 
for ip, and for some diagnostic tools. A smoothed version of the residuals is also proposed, leading 
to more appealing graphics. The main purpose of the residuals analysis is to check whether the 
fitted model is misspecified. 

As an illustration, let us assess the effect of a misspecified model to the tessellation simulated 
in top left of Figure[S] It actually corresponds to a sample from Model 3 where a = 0.05, B = +00, 
z = 100 and 9 ~ —0.5. But we will improperly fit a stationary Poisson process to this sample, 
then we will fit Model 2 (where the interaction relies on the Delaunay triangulation). Finally the 
correct Model 3 will be fitted for a sake of comparison. Figure [T^] represents the sample according 
to these three points of view. 
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Figure 13: Voronoi tessellation (left) and Delaunay tessellation (right) from the same point configuration 
(middle) , coming from a simulation of Model 3 (top right of Figure [B|) . 

We consider the simple case when -ip — 1. This corresponds to the so-called raw residuals, which 
have the following form in the stationary case: 
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To check the fitted model, we use the QQ-plot diagnostic presented in [1|. It consists in comparing 
the empirical quantiles of the fitted residuals to the empirical quantiles of bootstrapped residuals. 

If we fit a stationary Poisson process to the sample, we obtain an estimated intensity z = 833. 
The raw residuals, computed on squares A with side 0.01, are shown in top left of Figure [HI The 
same kind of residuals have then been computed on 100 simulated Poisson process with intensity 
z. A QQ-plot of these residuals with a 95%-confidence interval is shown on top right of Figure [T4j 
where the residuals of the original sample have been added (crosses) . An example of raw residuals 
from a simulated Poisson process is represented in the middle of this plot. It appears that the 
residuals of our sample do not behave as those from the simulated Poisson point processes. The 
stationary Poisson model is then misspecified. 

Similarly, if we fit Model 2 to the same sample, we obtain a — 0.055 and d = 4.49 when 
z — 1000. The raw residuals for this model are shown in bottom left of Figure [TH We have 
bootstrapped residuals from 100 simulated samples from Model 2 with the same parameters. One 
example of such residuals is shown on bottom middle. The QQ-plot, in bottom right, shows that 
the original sample does not seem to follow Model 2. 

Finally, Model 3 is fitted. The estimation gives a = 0.049, B = 97.8 and 9 = -0.56 when 
z = 100. The same plots as before are represented in Figure [151 According to the QQ-plot, one 
should not reject the fitted model. Let us remark that the estimation of 9 is rather bad for our 
sample: the error is —0.06 although other simulations shows that the standard deviation of the 
errors is about 0.02. This is the reason why the distribution of the residuals is on the edge of 
the confidence interval in the QQ-plot. The two residuals images represented on the left show 
that some residuals may be very negative on some squares (the black ones). This is confirmed by 
the dispersion of the lowest quantiles in the QQ-plot. Thus, the distribution of the residuals can 
not be Gaussian in this example. This differs from the asymptotic gaussianity of most residuals 
conjectured in [l|. 

AppendixA. 

AppendixA.l. Existence of Gibbs Delaunay-Voronoi tessellations. 

The existence results presented here are published in a more general setting in |[8i] and |9|] . They 
are slightly modified and simplified so that they suit better the setting of Gibbs Delaunay-Voronoi 
tessellations. We suppose that the energy functions have the forms ((H) or ([3]). The three following 
assumptions H1-H3 are sufficient to define the conditional densities /a in ^. 

HI Xo,(R2)^0. 

For every 7 in A^oo(R^) and every A in B{B?), 7 in 7W(R^) is called a (r, A)-modification of 7 
(with r > 0) if there exist distinct yi, 2/2, • ■ • , J/n in A satisfying |yi — ccij < r for every 1 < i < n 
(with 7A = {xi,X2,-..,Xn}) such that 7 = {yi, 2/2, •■-,?/«} U 7ac. 

H2 A^oo(K^) is a locally open set in A^(R^) which means that for every 7 in A^oo(R^), every A 
in B{S?) there exists rA(7) > such that any (rA(7), A)-modification of 7 is in A^oo(K^)- 

H3 The interactions Vi and V2 are stable which means that there exists a constant K > such 
that 

Vi > ~K and V2 > -K. 

Now let us give a collection of assumptions used in the proof of the existence of Gibbs Delaunay- 
Voronoi tessellations. For every i? > 0, we denote by 7^2 an infinite configuration in i?-equilateral 
position. That means that any triangle in 'Del{'yfi) is equilateral with length of sides equal to R 
(see the initial configuration presented in Figure [2| for an example). 

H4 There exist ATi > and K[ > such that, for all 7 £ Moo{^'^) and ah A e ^(R^), 

Card(7A) > i^iVo/(A) - K[. 
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Residuals when fitting Model 2 Simulated residuals from Model 2 QQplot from bootstrap 

Figure 14: Analysis of residuals for misspecified models. 
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Residuals when fitting Model 3 Simulated residuals from Model 3 QQplot from bootstrap 

Figure 15: Analysis of residuals for the correct model. 
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H5 There exist R > 0, < r < R/2 and A > such that for every A e B{M?), every (r,A)- 
inodification 7 of ^ji, and every T, T' € 'Del{'y) with T ^Dej T', 

Vi{T)<A and V2(r,T')<A; (A.l) 

Respectively, for every C, C" G Vor(7) with C ^vor C", 

Vi{C)<A and F2(C,C")<A (A.2) 

Now we are able to give a first existence theorem 

Theorem 1. There exists a stationary Gihbs Delaunay-Voronoi tessellation for any intensity 
v = z\ (z > Q) and any energy Junctions (£'A)AeB(R2) satisfying assumptions HI, H2, H3, H4 
and H5. 

Assumptions H4 and H5 can be substituted by the following one. 

H6 There exists A > Q such that for every r > we can find R > 2r such that for every 
A G B(R^) and every (r, A)-modification 7 of 77^ the property (jA.l|l or (|A.2p holds. 

We have the second following existence theorem. 

Theorem 2. There exists a stationary Gihhs Delaunay-Voronoi tessellation for any intensity 
v = zX (z > Q) and any energy functions (-Ea)agb(R2) satisfying the Assumptions HI, H2, H3 
and H6. 

The proofs of Theorems [T] and [5] can be found in [9] . They rely on entropy tools which are 
only available in the setting of stationary processes (i.e. v = zX). Concerning our three example 
models, the following corollary holds (the existence of Model 2 is also proved in [8|). 

Corollary 1. In the stationary case, i.e. when the intensity measure v is equal to zX, Gihbs 
Delaunay- Voronoi tessellations for models 1, 2 and 3 exist. 

Proof. First of all, assumptions HI, H2, H3 are obviously satisfied for the three models. Con- 
cerning Model 1, the existence is given by Theorem [51 Assumption H6 is proved by taking A = 
and R large enough with respect to r such that any (r, A)-modification 7 of ^b. have Delaunay 
triangles with angles larger than a. Concerning Models 2 and 3, the existence is given by Theorem 
[TJ Assumption H4 comes from the hardcore interaction which forces the cells to be not too large. 
The uniform bound in H5 is obvious if R and r are chosen such that any (r, A)-modification 7 of 
7i?isinXoo(K^). □ 

Appendix A. 2. Gonvergence of the algorithm 

The Birth-Death-Move algorithm used in this paper is presented in [22] page 115 where the 
convergence is proved in Proposition 7.7 if the associated Markov Chain is irreducible and aperiodic 
(see also [21|). In our setting, there is no problem with the aperiodicity since the probability that 
nothing happens during one step of the algorithm is positive. In general to prove the irreducibility, 
it is sufficient to point out that every configuration 7 is linked by a finite number of algorithm 
steps to the empty configuration. In our case, it is not possible because there is a strong hardcore 
interaction and so the connection with the empty configuration is in general false. So we need the 
connectivity assumption H7 based on the following definition. 

Definition 5. 7 and 7' in A^oo(K^) are connected if there exist n > and a sequence of configu- 
rations 7o, 71, ... , 7n-i, 7n "in A^oo(II^^) such that 70 = 7, 7n = 7' and, for each < i < n — 1, ji 
and 7i+i differ only by one step of the algorithm (a birth, a death or a move). 

H7 For every 7 and 7' in [0, 1]^ such that 7 and 7', defined in ([S]), are in 7V{oo(R^), then 7 and 
7' are connected. 
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The deterministic connectivity assumption H7 and the flexibihty assumption H2 on the space 
A^oo(lK^) ensure that for all configurations 7, 7' in A^oo(R^) and every r > 0, the algorithm may 
generate from 7, with a positive probability and a finite number of steps, a (r, [0, l]^)-modification 
of 7' (see the definition after HI). Irreducibility of the Markov chain follows and we have the 
following proposition. 

Proposition 3. Under the assumptions H2 and H7, the Birth- Death- Move algorithm presented 
in Section 3.2 converges to P (see Definitions^ in total variation norm for P-a.s. every initial 
condition. 

It seems difficult in our context to obtain rates of convergence, because the energy functions 
are not locally stable (the local stability means that \E\{j U {x}) — i?A(7)| is uniformly bounded 
with respect to A, 7 and x). Moreover the space A^oo(R^) may be very complicated since there is 
no upper bound in general for the number of steps n in assumption H7. 

Let us remark that assumption H7 is not easy to check. If H7 is not satisfied then the algorithm 
converges to the restriction of P on the connected component of the initial configuration 70 for 
the connection relation defined below (see definition [5]) . In this case, the limiting distribution may 
depend on the initial configuration. For Model 1, we can show that H7 is satisfied if a is small 
enough. We don't give the proof here but the scheme is essentially the same than in the following 
Proposition [4] which deals with Model 2. For Model 3, it is more complicated, we have not proved 
it but it seems satisfied if a and B are large enough. 

Proposition 4. In Model 2, if 2e < a < ^ then the algorithm presented in Sections 3.2 converges 
to P. 

Proof. According to Proposition [31 it suffices to show H2 and H7. Since H2 is obviously satisfied 
for Model 2, it remains to show H7. 

Let 7 and 7' be in [0, 1]^ such that 7 and 7' are in A^oo(IR^)- To simplify the notations, we 
say that 7 is in M.oo{^^) if 7 is in A^oc(R^)- We start the sequence by putting 70 = 7 and we 
construct the sequence 7^ by an algorithmic procedure. 

In a first step (called saturation) we add points until there does not exist any ball with radius 
£ without points. More precisely, we test if there exists x in [0, 1]^ such that 70 n B{x, e) = 0. If it 
is not the case, the saturation is finished. If it is the case we add the point a; by a birth-step action 
and we put 71 = 70 + x. Then, we test again if there exists x in [0, 1]^ such that 71 n B{x., e) — 0. 
If it is not the case the saturation is finished otherwise we put 72 = 71 + x. We go on like this 
until the saturation procedure stops which is always the case since [0,1]^ may contain only a finite 
number of points with a distance between them bigger than e. Let us remark that this construction 
produces configurations in A^oo(R^)- We denote by 7^,^ the saturated configuration of 7. 

In a second step, we add the points of 7' to 7^1 by the following way. Let a; be a point of 7'. 
By definition of the saturation, the configuration 7^1 + x is not in A^oo(R^) since there exists at 
least one point y in 7^^ such that |x — yj < £ {y is the periodic version of y such that y G B(x, e)). 
If this point y is unique we use a move-step action to move y to x. So we put 7mi+i — (7 ~ y) + 2;. 
If these points are non unique, they are removed (except one) by death-step actions and the last 
one is moved as above. We denote by 7^1+1, •■• ,7m2 this sequence and we remark that these 
configurations are in A^oo(R^) since 7^,^ is saturated and 2e < a < i. Now we saturate again the 
configuration 7,„2 as above and we add another point of 7' to 7^2 • We go on until we have added 
all the points of 7' and we denote by 7^3 the final configuration. 

It remains to remove the points of jma which are not in 7'. It is sufficient to apply several 
death-step actions since the obtained configurations are in 7V{oo(R^)- 

D 

Appendix A. 3. Consistency of the estimation procedure 

Let us suppose that the energy function, defined in (U} and ([S]), is parameterized by /3 and 9, 
and is denoted by E^' . We first need to distinguish properly the hardcore parameter /3 from the 
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other parameter 9. This is the purpose of the following assumption. 
SI: For all 7 G A4(M2), f^^. ^^^ ^ .^^^ ^j. j^jj ^ ^^j^j 5,/^ 

VA e 6(m2), £;^^''(7a,7aO < cx) ^^ i?^^''(7A,7A0 < ^- 

Under SI. the support of the energy is parameterized by /3 only, and not by 6, which confirms 
that /? is the hardcore parameter. This assumption is satisfied by Models 2 and 3 with /3 = (e, a) 
and P — (e, a, B) respectively. 

The strong consistency (under any stationary Gibbs measure) of /3 defined in pop and (z, 0) 
defined in (|13p are established in Theorem 2 in [10], under some regularity assumptions. These 
assumptions have been checked for Model 2 in Proposition 5 in [10]. Concerning Model 3, the 
assumptions could be checked in the same way excepted for assumption S3 involved in [lOi]. We 
have not succeeded to prove it but it seems true at least for B large enough. 
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